Aprendizado Supervisionado

O objetivo desse notebook é aplicar a Linguagem R em Ciência de Dados na resolução de problemas de classificação.

Classificando flores

A base de dados iris fornece medidas em centímetros das variáveis largura e comprimento da sépala e largura e comprimento da pétala, respectivamente, para 50 flores, de cada uma das 3 espécies de íris.

data('iris')

Podemos usar o scatterplot para entender o comportamento das classes

plot(iris[,1:4], col=iris$Species)

Devemos normalizar, uma vez que podemos utilizar algoritmos de classificação que são suscetíveis aos valores dos atributos.

iris = data.frame(scale(iris[,1:4]), Species=iris$Species)
summary(iris)

##   Sepal.Length       Sepal.Width       Petal.Length      Petal.Width     
##  Min.   :-1.86378   Min.   :-2.4258   Min.   :-1.5623   Min.   :-1.4422  
##  1st Qu.:-0.89767   1st Qu.:-0.5904   1st Qu.:-1.2225   1st Qu.:-1.1799  
##  Median :-0.05233   Median :-0.1315   Median : 0.3354   Median : 0.1321  
##  Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.67225   3rd Qu.: 0.5567   3rd Qu.: 0.7602   3rd Qu.: 0.7880  
##  Max.   : 2.48370   Max.   : 3.0805   Max.   : 1.7799   Max.   : 1.7064  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
## 

Depois de carregar a base. Precisamos particionar o conjunto de treinamento e teste. Para essa e outras bases vamos usar o holdout com 70% para treinamento e 30% para teste.

set.seed(123)
idx = sample(nrow(iris), 0.7*nrow(iris), replace=FALSE)
tran = iris[idx,]
test = iris[-idx,]

Podemos conferir se ocorreu contaminação.

intersect(rownames(tran), rownames(test))

## character(0)

Agora podemos induzir alguns modelos de Aprendizado de Máquina. Vamos começar pela Árvore de Decisão.

require('rpart')

## Loading required package: rpart

model = rpart(Species~., iris)

No caso da Árvore de Decisão, é possível plotar o modelo. Para isso vamos usar a biblioteca rattle.

library('rattle')

## Loading required package: tibble

## Loading required package: bitops

## Rattle: A free graphical interface for data science with R.
## Version 5.4.0 Copyright (c) 2006-2020 Togaware Pty Ltd.
## Type 'rattle()' to shake, rattle, and roll your data.

fancyRpartPlot(model, caption=NULL)

Também podemos predizer a classe das amostras de teste.

pred = predict(model, test, type='class')

E plotar o scatterplot para os dados de treinamento e teste.

par(mfrow=c(1,2))
plot(tran[,3:4], col=tran$Species, pch=as.numeric(tran$Species))
plot(test[,3:4], col=pred, pch=as.numeric(pred))

Podemos avaliar de forma mais objetiva o desempenho do modelo.

conf = table(test$Species, pred)
print(conf)

##             pred
##              setosa versicolor virginica
##   setosa         14          0         0
##   versicolor      0         18         0
##   virginica       0          1        12

Inclusive, calcular a taxa de acerto.

acc = sum(diag(conf))/sum(conf)
print(acc)

## [1] 0.9777778

Existem outras medidas de avaliação. Recomendo olhar os pacotes ‘caret’ e ‘mlr’.

Agora podemos induzir outros modelos de Aprendizado de Máquina. Para isso precisamos carregar algumas bibliotecas.

require('kknn')

## Loading required package: kknn

require('e1071')

## Loading required package: e1071

require('randomForest')

## Loading required package: randomForest

## randomForest 4.6-14

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:rattle':
## 
##     importance

Gerando os modelos de Aprendizado de Máquina.

model1 = kknn(Species~., tran, test, k=3)
model2 = svm(Species~., tran, kernel="radial")
model3 = randomForest(Species~., tran)

Calculando a predição.

pred1 = model1$fitted.values
pred2 = predict(model2, test, type='class')
pred3 = predict(model3, test)

Calculando o desempenho baseado na acurácia.

acuracia <- function(orig, pred) {
  acc = table(orig, pred)
  sum(diag(acc))/sum(acc)  
} 

cat('A acurácia do 3NN é', acuracia(test$Species, pred1))

## A acurácia do 3NN é 0.9333333

cat('A acurácia do SVM é', acuracia(test$Species, pred2))

## A acurácia do SVM é 0.9777778

cat('A acurácia do RF é', acuracia(test$Species, pred3))

## A acurácia do RF é 0.9777778

No entanto, o holdout não é a melhor alternativa. A validação cruzada estratificada é mais interessante. Podemos ver na imagem um exemplo de seu funcionamento.

k-fold cross validation (Adaptado da Wikepedia)

Podemos implementar :)

kfolds <- function(y, folds) {

  names(y) =  1:length(y)
  index = lapply(1:nlevels(y), function(i) {
    rep(1:folds, length.out=length(y[y == levels(y)[i]]))
  })

  index = unlist(index)
  folds = lapply(1:folds, function(i) {
    as.numeric(names(y[index == i]))
  })

  return(folds)
}

E agora podemos usar.

idx = kfolds(iris$Species, 10)

acc = lapply(idx, function(i){
  model = randomForest(Species~., iris[-i,])
  pred = predict(model3, iris[i,])
  acuracia(iris[i,]$Species, pred)
})

cat(mean(unlist(acc)), '+-', sd(unlist(acc)))

## 0.9933333 +- 0.02108185

Classificando dígitos

Um experimento interessante é a tarefa de classificar dígitos utilizando o algoritmo de Redes Neurais Artificiais. Para isso precisamos relizar os experimentos em duas fases: (1) treinar a rede para reconhecer os dígitos e (2) utilizá-la para classificar novos dígitos. Para isso precisamos carregar as blbiotecas ‘tensorflow’ e ‘keras’.

require('tensorflow')

## Loading required package: tensorflow

require('keras')

## Loading required package: keras

A base de dados de dígitos se chama MNIST e pode ser carregada com a função dataset_mnist. Ela contém 60000 amostras dos dígitos 0 até 9. Cada dígito é uma imagem de 28 x 28 pixels.

mnist = dataset_mnist()

A função show_digit recebe um dígito na forma de uma matriz 28 x 28 e imprime o dígito na forma de imagem.

show_digit <- function(numero) {
  image(numero[1:28,28:1])
}

Podemos utilizar a função show_digit para plotar uma amostra de todos os dígitos presentes na base de dados.

par(mfrow=c(2,5))
plot = lapply(0:9, function(i) {
  aux = mnist$train$x[mnist$train$y == i,,][1,,]
  show_digit(t(aux))
})

Normalizando os dados de treinamento e teste.

mnist$train$x = mnist$train$x/255
mnist$test$x = mnist$test$x/255

Construindo uma Rede Neural Artificial usando Keras. Essa rede é composta de uma camanda de entrada flatten (achatada de 28 x 28 neurônios), uma camada oculta dense (128 neurônios com função de ativação linear) e uma camada de saída dense (10 neurônios com função de ativação exponencial normalizada).

modelo = keras_model_sequential() %>% 
  layer_flatten(input_shape = c(28, 28)) %>% 
  layer_dense(units = 128, activation = "relu") %>% 
  layer_dense(10, activation = "softmax")

Informações sobre a Rede Neural Artificial.

summary(modelo)

## Model: "sequential"
## ________________________________________________________________________________
## Layer (type)                        Output Shape                    Param #     
## ================================================================================
## flatten (Flatten)                   (None, 784)                     0           
## ________________________________________________________________________________
## dense_1 (Dense)                     (None, 128)                     100480      
## ________________________________________________________________________________
## dense (Dense)                       (None, 10)                      1290        
## ================================================================================
## Total params: 101,770
## Trainable params: 101,770
## Non-trainable params: 0
## ________________________________________________________________________________

Podemos representar graficamente:

library(deepviz)
library(magrittr)
modelo %>% plot_model()

## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

Ainda precisamos definir o algoritmo de treinamento, a função de avaliação e a medida de desempenho que vamos utilizar. Podemos fazer isso com a função compile.

modelo %>% 
  compile(
    loss = "sparse_categorical_crossentropy",
    optimizer = "adam",
    metrics = "accuracy"
  )

Agora podemos treinar a rede.

historico = modelo %>% 
  fit(
    x = mnist$train$x, y = mnist$train$y,
    epochs = 10,
    validation_split = 0.3,
    verbose = 2
  )

É interessante medir o erro e o desempenho da rede ao longo das épocas no conjunto de treinamento e validação.

plot(historico)

## `geom_smooth()` using formula 'y ~ x'

Agora que sabemos que a rede aprendeu a reconhecer dígitos, podemos verificar seu desempenho no conjunto de teste. Para isso precisamos aplicar essas amostras nunca antes vistas na entrada da rede e obter os valores de saída.

predicao = predict(modelo, mnist$test$x)
head(predicao)

##              [,1]         [,2]         [,3]         [,4]         [,5]
## [1,] 3.084949e-07 5.754219e-12 4.037815e-06 5.538938e-04 5.470358e-13
## [2,] 3.502256e-10 1.211353e-06 9.999982e-01 1.676618e-07 3.619702e-20
## [3,] 1.660483e-06 9.990133e-01 2.579817e-04 9.216476e-06 1.853301e-05
## [4,] 9.999995e-01 3.414662e-11 4.631823e-08 1.901913e-11 3.870555e-13
## [5,] 3.807434e-09 1.649097e-10 3.516960e-08 8.413788e-11 9.989899e-01
## [6,] 2.665162e-08 9.998041e-01 2.826536e-07 1.621908e-07 1.830018e-06
##              [,6]         [,7]         [,8]         [,9]        [,10]
## [1,] 9.990508e-09 4.179083e-16 9.994093e-01 3.228376e-06 2.929628e-05
## [2,] 2.507383e-07 1.856738e-09 7.926855e-17 6.413249e-08 1.516677e-14
## [3,] 1.012086e-06 1.325644e-05 1.331137e-04 5.498987e-04 2.033516e-06
## [4,] 1.341752e-09 1.191723e-09 1.357380e-07 2.338699e-13 3.172412e-07
## [5,] 2.498284e-11 8.186302e-09 8.766361e-07 8.853513e-09 1.009312e-03
## [6,] 7.383683e-10 9.211697e-09 1.890462e-04 4.439606e-06 1.114614e-07

Cada neurônio de saída da rede irá retornar um valor. Aquele neurônio com maior valor deve indicar a classe com maior proababilidade.

predicao = apply(predicao, 1, which.max)
predicao[1:6]

## [1] 8 3 2 1 5 2

Calculando o desempenho da rede no conjunto de teste.

acc = acuracia(mnist$test$y, predicao)
cat('A acurácia do modelo no conjunto de teste eh:', acc)

## A acurácia do modelo no conjunto de teste eh: 0.9757

Agora podemos usar os modelos de Árvore de Decisão, kNN, SVM e Random Forest para comparar com o desempenho com a Rede Neural Artificial. Primeiro vamos construir uma base tabular com a mnist.

tran = data.frame(mnist$train$x, y=as.factor(mnist$train$y))
test = data.frame(mnist$test$x, y=as.factor(mnist$test$y))

Agora podemos induzir outros modelos.

model = rpart(y~., tran)
pred = predict(model, test, type='class')

cat('A acurácia da ANN é:', acc)

## A acurácia da ANN é: 0.9757

cat('A acurácia da DT é', acuracia(test$y, pred))

## A acurácia da DT é 0.6196